Raincloud plots of benchmark and simulation data under extremes

Overview

This notebook will guide you through the code to plot raincloud plots of the data under climatic extremes in R. The following figures will be computed:

• Figure 3: raincloud plots for relative detrended yield values for aggregated crops

• Supplementary figures 10-13: raincloud plots for relative detrended yield values for each crop (mai, soy, wwh, ri1)

• Supplementary figures 14-18: raincloud plots for relative detrended yield values filtered for top 5 producing countries

• Supplementary figures 19-23: raincloud plots for absolute detrended yield values

• Supplementary figures 23-27: raincloud plots for model errors using relative detrended yield values

• Supplementary figures 28-32: raincloud plots for model errors using absolute detrended yield values

The fully processed and filtered data files, available under GGCMI-validation/data/processed/figure_ready_data are required to produce the figures.

We load the necessary libraries

Code

library(tidyverse)   # Includes dplyr, ggplot2, readr, etc.
library(sjPlot)  # to set and configure plotting theme
library(rnaturalearth) # to create base world map with countries

We configure the plotting theme based on the data visualization blog from Cédric Scherer: https://www.cedricscherer.com/

Code

theme_Scherer <- function (base_size = 12, base_family = "Helvetica") {
  half_line <- base_size/2
  theme(
    line = element_line(color = "black", linewidth = .5,
                        linetype = 1, lineend = "butt"),
    rect = element_rect(fill = "white", color = "black",
                        linewidth = .5, linetype = 1),
    text = element_text(family = base_family, face = "plain",
                        color = "black", size = base_size,
                        lineheight = .9, hjust = .5, vjust = .5,
                        angle = 0, margin = margin(), debug = FALSE),
    axis.line = element_blank(),
    axis.line.x = NULL,
    axis.line.y = NULL,
    axis.text = element_text(size = base_size * 1.1, color = "gray30"),
    axis.text.x = element_text(margin = margin(t = .8 * half_line/2),
                               vjust = 1),
    axis.text.x.top = element_text(margin = margin(b = .8 * half_line/2),
                                   vjust = 0),
    axis.text.y = element_text(margin = margin(r = .8 * half_line/2),
                               hjust = 1),
    axis.text.y.right = element_text(margin = margin(l = .8 * half_line/2),
                                     hjust = 0),
    axis.ticks = element_line(color = "gray30", linewidth = .7),
    axis.ticks.length = unit(half_line / 1.5, "pt"),
    axis.ticks.length.x = NULL,
    axis.ticks.length.x.top = NULL,
    axis.ticks.length.x.bottom = NULL,
    axis.ticks.length.y = NULL,
    axis.ticks.length.y.left = NULL,
    axis.ticks.length.y.right = NULL,
    axis.title.x = element_text(margin = margin(t = half_line),
                                vjust = 1, size = base_size * 1.3,
                                face = "bold"),
    axis.title.x.top = element_text(margin = margin(b = half_line),
                                    vjust = 0),
    axis.title.y = element_text(angle = 90, vjust = 1,
                                margin = margin(r = half_line),
                                size = base_size * 1.3, face = "bold"),
    axis.title.y.right = element_text(angle = -90, vjust = 0,
                                      margin = margin(l = half_line)),
    legend.background = element_rect(color = NA),
    legend.spacing = unit(.4, "cm"),
    legend.spacing.x = NULL,
    legend.spacing.y = NULL,
    legend.margin = margin(.2, .2, .2, .2, "cm"),
    legend.key = element_rect(fill = "gray95", color = "white"),
    legend.key.size = unit(1.2, "lines"),
    legend.key.height = NULL,
    legend.key.width = NULL,
    legend.text = element_text(size = rel(.8)),
    legend.text.align = NULL,
    legend.title = element_text(hjust = 0),
    legend.title.align = NULL,
    legend.position = "right",
    legend.direction = NULL,
    legend.justification = "center",
    legend.box = NULL,
    legend.box.margin = margin(0, 0, 0, 0, "cm"),
    legend.box.background = element_blank(),
    legend.box.spacing = unit(.4, "cm"),
    panel.background = element_rect(fill = "white", color = NA),
    panel.border = element_rect(color = "gray30",
                                fill = NA, linewidth = .7),
    panel.grid.major = element_line(color = "gray90", linewidth = 1),
    panel.grid.minor = element_line(color = "gray90", linewidth = .5,
                                    linetype = "dashed"),
    panel.spacing = unit(base_size, "pt"),
    panel.spacing.x = NULL,
    panel.spacing.y = NULL,
    panel.ontop = FALSE,
    strip.background = element_rect(fill = "white", color = "gray30"),
    strip.text = element_text(color = "black", size = base_size),
    strip.text.x = element_text(margin = margin(t = half_line,
                                                b = half_line)),
    strip.text.y = element_text(angle = -90,
                                margin = margin(l = half_line,
                                                r = half_line)),
    strip.text.y.left = element_text(angle = 90),
    strip.placement = "inside",
    strip.placement.x = NULL,
    strip.placement.y = NULL,
    strip.switch.pad.grid = unit(0.1, "cm"),
    strip.switch.pad.wrap = unit(0.1, "cm"),
    plot.background = element_rect(color = NA),
    plot.title = element_text(size = base_size * 1.8, hjust = .5,
                              vjust = 1, face = "bold",
                              margin = margin(b = half_line * 1.2)),
    plot.title.position = "panel",
    plot.subtitle = element_text(size = base_size * 1.3,
                                 hjust = .5, vjust = 1,
                                 margin = margin(b = half_line * .9)),
    plot.caption = element_text(size = rel(0.9), hjust = 1, vjust = 1,
                                margin = margin(t = half_line * .9)),
    plot.caption.position = "panel",
    plot.tag = element_text(size = rel(1.2), hjust = .5, vjust = .5),
    plot.tag.position = "topleft",
    plot.margin = margin(rep(base_size, 4)),
    complete = TRUE
  )
}
set_theme(theme_Scherer())

Then we can load our data:

Code

base_path <- "C:/Users/kobed/Documents/TOAD-PIKlegacy" # Where is the code repo stored?

# List of crops
crops <- c("mai", "wwh", "ri1", "soy", "aggr")

# Read each .rds file into a named list
crop_data <- lapply(crops, function(crop) {
  readRDS(file.path(base_path, paste0("extremes_", crop, ".rds")))
})

names(crop_data) <- crops

Figure 3 – rainclouds for aggregated crops

First we prepare the data

Code

dat_long <- crop_data[["aggr"]] %>% 
  dplyr::select(lat, lon, year, divtrend_sim, divtrend_obs, model, climate_extreme) %>% 
mutate(divtrend_sim_perc = (divtrend_sim - 1)*100, # Transform to percentages
       divtrend_obs_perc = (divtrend_obs - 1)*100
      ) %>%
  pivot_longer(
    cols = c(divtrend_obs_perc, divtrend_sim_perc),  # Adjust to the actual column names you want to pivot
    names_to = "type",
    values_to = "divtrend"
  ) %>% 
  mutate(
    model = case_when(
      type == "divtrend_obs_perc" ~ "benchmark",  
      type == "divtrend_sim_perc" ~ model
    )
  ) %>% 
  distinct() 

Then we can calculate the proportion of underestimation and the median detrended yield for the benchmark data

Code

underestimation_rate <- crop_data[["aggr"]] %>%
  mutate(error_perc = (divtrend_obs - 1)*100 - (divtrend_sim - 1)*100) %>%
  mutate(error_perc = ifelse(is.infinite(error_perc), -1000, error_perc)) %>%  # Replace Inf with -1000
  group_by(model, climate_extreme) %>%
  summarise(under_rate = mean(error_perc < 0) * 100)  # Calculate % of values below 0

benchmark_medians <- dat_long %>%
  filter(model == "benchmark") %>%
  group_by(climate_extreme) %>%
  summarise(median_yield = median(divtrend, na.rm = TRUE))

Prepare the plot by selecting the name for the benchmark and creating a fixed order to list the models

Code

# Define the name of the benchmark distribution
benchmark_name <- "benchmark"  

# Model order
order <- rev(c(
  "benchmark",
  "ensemble",
  "acea", 
  "simplace-lintul5",
  "dssat-pythia",
  "pepic",
  "cygma1p74",
   "ldndc",
  "pdssat",
  "epic-iiasa", 
  "lpj-guess", 
  "promet",
  "lpjml",
  "isam", 
  "crover"
))

Now we can create the raincloud plot

Code

# Create the raincloud plot
raincloud_perc <- dat_long %>%
    mutate(model = factor(model, levels = order)) %>%  # Order the models using the defined order
   filter(is.finite(divtrend), abs(divtrend) < 1000) %>%  # Filter out Inf values and keep within range
  ggplot(aes(x = divtrend, y = model, fill = model == benchmark_name)) +
  
  
  # Add half-eye plot
  ggdist::stat_halfeye(aes(fill = model == benchmark_name), adjust = 0.7, 
                       width = 0.7,               # Increase distribution width
                       justification = -0.1,    # Move half-eye plot further from boxplot
                       .width = c(0.3, 0), 
                       point_colour = "black", 
                       alpha = 0.7) +
  
    
  scale_fill_manual(values = c("TRUE" = "#00A087FF",  
                             "FALSE" = "#7E6148FF"))+ 
  
  # Add boxplots (smaller width to give more space to distribution)
  geom_boxplot(width = 0.12, outlier.shape = NA, alpha = 0.7, 
               fill = "white", color = "black", size = 1) +
  
   # Add black vertical line for the benchmark median
  geom_vline(data = benchmark_medians, aes(xintercept = median_yield), 
              color = "black", size = 2) +
 
  # Add text annotations for the underestimation rate (place them to the right of the plot)
  geom_text(data = underestimation_rate, 
           aes(x = 40, y = model, label = paste0(round(under_rate, 1), "%")),  # Round percentage to 1 decimal
            vjust = -2, size = 5, color = "black") +  # Adjust text size and positioning
  
  # Labels and annotations
  labs(title = "", 
       x = "Yield anomalies in %", 
       y = "") +
  
  # Add facets for different climate extremes
  facet_wrap(~ climate_extreme, scales = "free_x", ncol = 2) +

  
  # Custom x-axis breaks
  scale_x_continuous(breaks = seq(-100, 100, by = 20)) +
  
   # Truncate the x-axis to focus on the range
  coord_cartesian(xlim = c(-100, 100)) +
  
   theme(legend.position = "none", 
        axis.title.x = element_text(size = 18),          # Increase x-axis title size
        axis.text.x = element_text(size = 18),           # Increase x-axis text size
        axis.text.y = element_text(size = 18),           # Increase y-axis text size
        strip.text = element_text(size = 18))            # Increase facet label size

raincloud_perc

Code

# Save the figure 
ggsave(file.path(base_path, "fig4_RaincloudExtremes.png"), raincloud_perc, width = 18, height = 20, dpi = 300)

Supplementary figures 10-13 – rainclouds per crop

RICE

Code

dat_long <- crop_data[["ri1"]] %>% 
  dplyr::select(lat, lon, year, divtrend_sim, divtrend_obs, model, climate_extreme) %>% 
mutate(divtrend_sim_perc = (divtrend_sim - 1)*100, # Transform to percentages
       divtrend_obs_perc = (divtrend_obs - 1)*100
      ) %>%
  pivot_longer(
    cols = c(divtrend_obs_perc, divtrend_sim_perc),  # Adjust to the actual column names you want to pivot
    names_to = "type",
    values_to = "divtrend"
  ) %>% 
  mutate(
    model = case_when(
      type == "divtrend_obs_perc" ~ "benchmark",  
      type == "divtrend_sim_perc" ~ model
    )
  ) %>% 
  distinct() 

underestimation_rate <- crop_data[["ri1"]] %>%
  mutate(error_perc = (divtrend_obs - 1)*100 - (divtrend_sim - 1)*100) %>%
  mutate(error_perc = ifelse(is.infinite(error_perc), -1000, error_perc)) %>%  # Replace Inf with -1000
  group_by(model, climate_extreme) %>%
  summarise(under_rate = mean(error_perc < 0) * 100)  # Calculate % of values below 0

benchmark_medians <- dat_long %>%
  filter(model == "benchmark") %>%
  group_by(climate_extreme) %>%
  summarise(median_yield = median(divtrend, na.rm = TRUE))

# Create the raincloud plot
raincloud_perc <- dat_long %>%
    mutate(model = factor(model, levels = order)) %>%  # Order the models using the defined order
   filter(is.finite(divtrend), abs(divtrend) < 1000) %>%  # Filter out Inf values and keep within range
  ggplot(aes(x = divtrend, y = model, fill = model == benchmark_name)) +
  
  
  # Add half-eye plot
  ggdist::stat_halfeye(aes(fill = model == benchmark_name), adjust = 0.7, 
                       width = 0.7,               # Increase distribution width
                       justification = -0.1,    # Move half-eye plot further from boxplot
                       .width = c(0.3, 0), 
                       point_colour = "black", 
                       alpha = 0.7) +
  
    
  scale_fill_manual(values = c("TRUE" = "#00A087FF",  
                             "FALSE" = "#7E6148FF"))+ 
  
  # Add boxplots (smaller width to give more space to distribution)
  geom_boxplot(width = 0.12, outlier.shape = NA, alpha = 0.7, 
               fill = "white", color = "black", size = 1) +
  
   # Add black vertical line for the benchmark median
  geom_vline(data = benchmark_medians, aes(xintercept = median_yield), 
              color = "black", size = 2) +
 
  # Add text annotations for the underestimation rate (place them to the right of the plot)
  geom_text(data = underestimation_rate, 
           aes(x = 40, y = model, label = paste0(round(under_rate, 1), "%")),  # Round percentage to 1 decimal
            vjust = -2, size = 5, color = "black") +  # Adjust text size and positioning
  
  # Labels and annotations
  labs(title = "", 
       x = "Yield anomalies in %", 
       y = "") +
  
  # Add facets for different climate extremes
  facet_wrap(~ climate_extreme, scales = "free_x", ncol = 2) +

  
  # Custom x-axis breaks
  scale_x_continuous(breaks = seq(-100, 100, by = 20)) +
  
   # Truncate the x-axis to focus on the range
  coord_cartesian(xlim = c(-100, 100)) +
  
   theme(legend.position = "none", 
        axis.title.x = element_text(size = 18),          # Increase x-axis title size
        axis.text.x = element_text(size = 18),           # Increase x-axis text size
        axis.text.y = element_text(size = 18),           # Increase y-axis text size
        strip.text = element_text(size = 18))            # Increase facet label size

raincloud_perc

Code

# Save the figure 
ggsave(file.path(base_path, "suppfig10_RaincloudExtremes.png"), raincloud_perc, width = 18, height = 20, dpi = 300)

WHEAT

Code

dat_long <- crop_data[["wwh"]] %>% 
  dplyr::select(lat, lon, year, divtrend_sim, divtrend_obs, model, climate_extreme) %>% 
mutate(divtrend_sim_perc = (divtrend_sim - 1)*100, # Transform to percentages
       divtrend_obs_perc = (divtrend_obs - 1)*100
      ) %>%
  pivot_longer(
    cols = c(divtrend_obs_perc, divtrend_sim_perc),  # Adjust to the actual column names you want to pivot
    names_to = "type",
    values_to = "divtrend"
  ) %>% 
  mutate(
    model = case_when(
      type == "divtrend_obs_perc" ~ "benchmark",  
      type == "divtrend_sim_perc" ~ model
    )
  ) %>% 
  distinct() 

underestimation_rate <- crop_data[["wwh"]] %>%
  mutate(error_perc = (divtrend_obs - 1)*100 - (divtrend_sim - 1)*100) %>%
  mutate(error_perc = ifelse(is.infinite(error_perc), -1000, error_perc)) %>%  # Replace Inf with -1000
  group_by(model, climate_extreme) %>%
  summarise(under_rate = mean(error_perc < 0) * 100)  # Calculate % of values below 0

benchmark_medians <- dat_long %>%
  filter(model == "benchmark") %>%
  group_by(climate_extreme) %>%
  summarise(median_yield = median(divtrend, na.rm = TRUE))

# Create the raincloud plot
raincloud_perc <- dat_long %>%
    mutate(model = factor(model, levels = order)) %>%  # Order the models using the defined order
   filter(is.finite(divtrend), abs(divtrend) < 1000) %>%  # Filter out Inf values and keep within range
  ggplot(aes(x = divtrend, y = model, fill = model == benchmark_name)) +
  
  
  # Add half-eye plot
  ggdist::stat_halfeye(aes(fill = model == benchmark_name), adjust = 0.7, 
                       width = 0.7,               # Increase distribution width
                       justification = -0.1,    # Move half-eye plot further from boxplot
                       .width = c(0.3, 0), 
                       point_colour = "black", 
                       alpha = 0.7) +
  
    
  scale_fill_manual(values = c("TRUE" = "#00A087FF",  
                             "FALSE" = "#7E6148FF"))+ 
  
  # Add boxplots (smaller width to give more space to distribution)
  geom_boxplot(width = 0.12, outlier.shape = NA, alpha = 0.7, 
               fill = "white", color = "black", size = 1) +
  
   # Add black vertical line for the benchmark median
  geom_vline(data = benchmark_medians, aes(xintercept = median_yield), 
              color = "black", size = 2) +
 
  # Add text annotations for the underestimation rate (place them to the right of the plot)
  geom_text(data = underestimation_rate, 
           aes(x = 40, y = model, label = paste0(round(under_rate, 1), "%")),  # Round percentage to 1 decimal
            vjust = -2, size = 5, color = "black") +  # Adjust text size and positioning
  
  # Labels and annotations
  labs(title = "", 
       x = "Yield anomalies in %", 
       y = "") +
  
  # Add facets for different climate extremes
  facet_wrap(~ climate_extreme, scales = "free_x", ncol = 2) +

  
  # Custom x-axis breaks
  scale_x_continuous(breaks = seq(-100, 100, by = 20)) +
  
   # Truncate the x-axis to focus on the range
  coord_cartesian(xlim = c(-100, 100)) +
  
   theme(legend.position = "none", 
        axis.title.x = element_text(size = 18),          # Increase x-axis title size
        axis.text.x = element_text(size = 18),           # Increase x-axis text size
        axis.text.y = element_text(size = 18),           # Increase y-axis text size
        strip.text = element_text(size = 18))            # Increase facet label size

raincloud_perc

Code

# Save the figure 
ggsave(file.path(base_path, "suppfig11_RaincloudExtremes.png"), raincloud_perc, width = 18, height = 20, dpi = 300)

SOY

Code

dat_long <- crop_data[["soy"]] %>% 
  dplyr::select(lat, lon, year, divtrend_sim, divtrend_obs, model, climate_extreme) %>% 
mutate(divtrend_sim_perc = (divtrend_sim - 1)*100, # Transform to percentages
       divtrend_obs_perc = (divtrend_obs - 1)*100
      ) %>%
  pivot_longer(
    cols = c(divtrend_obs_perc, divtrend_sim_perc),  # Adjust to the actual column names you want to pivot
    names_to = "type",
    values_to = "divtrend"
  ) %>% 
  mutate(
    model = case_when(
      type == "divtrend_obs_perc" ~ "benchmark",  
      type == "divtrend_sim_perc" ~ model
    )
  ) %>% 
  distinct() 

underestimation_rate <- crop_data[["soy"]] %>%
  mutate(error_perc = (divtrend_obs - 1)*100 - (divtrend_sim - 1)*100) %>%
  mutate(error_perc = ifelse(is.infinite(error_perc), -1000, error_perc)) %>%  # Replace Inf with -1000
  group_by(model, climate_extreme) %>%
  summarise(under_rate = mean(error_perc < 0) * 100)  # Calculate % of values below 0

benchmark_medians <- dat_long %>%
  filter(model == "benchmark") %>%
  group_by(climate_extreme) %>%
  summarise(median_yield = median(divtrend, na.rm = TRUE))

# Create the raincloud plot
raincloud_perc <- dat_long %>%
    mutate(model = factor(model, levels = order)) %>%  # Order the models using the defined order
   filter(is.finite(divtrend), abs(divtrend) < 1000) %>%  # Filter out Inf values and keep within range
  ggplot(aes(x = divtrend, y = model, fill = model == benchmark_name)) +
  
  
  # Add half-eye plot
  ggdist::stat_halfeye(aes(fill = model == benchmark_name), adjust = 0.7, 
                       width = 0.7,               # Increase distribution width
                       justification = -0.1,    # Move half-eye plot further from boxplot
                       .width = c(0.3, 0), 
                       point_colour = "black", 
                       alpha = 0.7) +
  
    
  scale_fill_manual(values = c("TRUE" = "#00A087FF",  
                             "FALSE" = "#7E6148FF"))+ 
  
  # Add boxplots (smaller width to give more space to distribution)
  geom_boxplot(width = 0.12, outlier.shape = NA, alpha = 0.7, 
               fill = "white", color = "black", size = 1) +
  
   # Add black vertical line for the benchmark median
  geom_vline(data = benchmark_medians, aes(xintercept = median_yield), 
              color = "black", size = 2) +
 
  # Add text annotations for the underestimation rate (place them to the right of the plot)
  geom_text(data = underestimation_rate, 
           aes(x = 40, y = model, label = paste0(round(under_rate, 1), "%")),  # Round percentage to 1 decimal
            vjust = -2, size = 5, color = "black") +  # Adjust text size and positioning
  
  # Labels and annotations
  labs(title = "", 
       x = "Yield anomalies in %", 
       y = "") +
  
  # Add facets for different climate extremes
  facet_wrap(~ climate_extreme, scales = "free_x", ncol = 2) +

  
  # Custom x-axis breaks
  scale_x_continuous(breaks = seq(-100, 100, by = 20)) +
  
   # Truncate the x-axis to focus on the range
  coord_cartesian(xlim = c(-100, 100)) +
  
   theme(legend.position = "none", 
        axis.title.x = element_text(size = 18),          # Increase x-axis title size
        axis.text.x = element_text(size = 18),           # Increase x-axis text size
        axis.text.y = element_text(size = 18),           # Increase y-axis text size
        strip.text = element_text(size = 18))            # Increase facet label size

raincloud_perc

Code

# Save the figure 
ggsave(file.path(base_path, "suppfig12_RaincloudExtremes.png"), raincloud_perc, width = 18, height = 20, dpi = 300)

MAIZE

Code

dat_long <- crop_data[["mai"]] %>% 
  dplyr::select(lat, lon, year, divtrend_sim, divtrend_obs, model, climate_extreme) %>% 
mutate(divtrend_sim_perc = (divtrend_sim - 1)*100, # Transform to percentages
       divtrend_obs_perc = (divtrend_obs - 1)*100
      ) %>%
  pivot_longer(
    cols = c(divtrend_obs_perc, divtrend_sim_perc),  # Adjust to the actual column names you want to pivot
    names_to = "type",
    values_to = "divtrend"
  ) %>% 
  mutate(
    model = case_when(
      type == "divtrend_obs_perc" ~ "benchmark",  
      type == "divtrend_sim_perc" ~ model
    )
  ) %>% 
  distinct() 

underestimation_rate <- crop_data[["mai"]] %>%
  mutate(error_perc = (divtrend_obs - 1)*100 - (divtrend_sim - 1)*100) %>%
  mutate(error_perc = ifelse(is.infinite(error_perc), -1000, error_perc)) %>%  # Replace Inf with -1000
  group_by(model, climate_extreme) %>%
  summarise(under_rate = mean(error_perc < 0) * 100)  # Calculate % of values below 0

benchmark_medians <- dat_long %>%
  filter(model == "benchmark") %>%
  group_by(climate_extreme) %>%
  summarise(median_yield = median(divtrend, na.rm = TRUE))

# Create the raincloud plot
raincloud_perc <- dat_long %>%
    mutate(model = factor(model, levels = order)) %>%  # Order the models using the defined order
   filter(is.finite(divtrend), abs(divtrend) < 1000) %>%  # Filter out Inf values and keep within range
  ggplot(aes(x = divtrend, y = model, fill = model == benchmark_name)) +
  
  
  # Add half-eye plot
  ggdist::stat_halfeye(aes(fill = model == benchmark_name), adjust = 0.7, 
                       width = 0.7,               # Increase distribution width
                       justification = -0.1,    # Move half-eye plot further from boxplot
                       .width = c(0.3, 0), 
                       point_colour = "black", 
                       alpha = 0.7) +
  
    
  scale_fill_manual(values = c("TRUE" = "#00A087FF",  
                             "FALSE" = "#7E6148FF"))+ 
  
  # Add boxplots (smaller width to give more space to distribution)
  geom_boxplot(width = 0.12, outlier.shape = NA, alpha = 0.7, 
               fill = "white", color = "black", size = 1) +
  
   # Add black vertical line for the benchmark median
  geom_vline(data = benchmark_medians, aes(xintercept = median_yield), 
              color = "black", size = 2) +
 
  # Add text annotations for the underestimation rate (place them to the right of the plot)
  geom_text(data = underestimation_rate, 
           aes(x = 40, y = model, label = paste0(round(under_rate, 1), "%")),  # Round percentage to 1 decimal
            vjust = -2, size = 5, color = "black") +  # Adjust text size and positioning
  
  # Labels and annotations
  labs(title = "", 
       x = "Yield anomalies in %", 
       y = "") +
  
  # Add facets for different climate extremes
  facet_wrap(~ climate_extreme, scales = "free_x", ncol = 2) +

  
  # Custom x-axis breaks
  scale_x_continuous(breaks = seq(-100, 100, by = 20)) +
  
   # Truncate the x-axis to focus on the range
  coord_cartesian(xlim = c(-100, 100)) +
  
   theme(legend.position = "none", 
        axis.title.x = element_text(size = 18),          # Increase x-axis title size
        axis.text.x = element_text(size = 18),           # Increase x-axis text size
        axis.text.y = element_text(size = 18),           # Increase y-axis text size
        strip.text = element_text(size = 18))            # Increase facet label size

raincloud_perc

Code

# Save the figure 
ggsave(file.path(base_path, "suppfig13_RaincloudExtremes.png"), raincloud_perc, width = 18, height = 20, dpi = 300)

Supplementary figures 14-18 – rainclouds for top producers

We need the unfiltered benchmark data to find out who the biggest producers are for this dataset.

Code

# Read each .rds file into a named list
crop_datageneral <- lapply(crops, function(crop) {
  readRDS(file.path(base_path, paste0("general_", crop, ".rds")))
})

general_sim_aggr <- readRDS(file.path(base_path, "general_sim_aggr.rds"))

names(crop_datageneral) <- crops

AGGREGATED

First we calculate who the 5 largest producing countries are and also calculate how much tons they produce and what the global share is (percentage)

Code

crop_datageneral[["aggr"]] %>% 
  left_join(general_sim_aggr, by = c("lon", "lat", "year")) %>%
  mutate(prod = yield * total_area) %>%
  group_by(ctr) %>%
  summarise(ctr_prod = sum(prod, na.rm = TRUE), .groups = "drop") %>%
  mutate(global_prod = sum(ctr_prod),
         global_share = (ctr_prod / global_prod) *100) %>%
  slice_max(order_by = ctr_prod, n = 5) %>%
  dplyr::select(ctr, ctr_prod, global_share)

# A tibble: 5 × 3
  ctr   ctr_prod global_share
  <chr>    <dbl>        <dbl>
1 CHN    5.63e15        23.5 
2 USA    5.53e15        23.1 
3 BRA    1.71e15         7.14
4 ARG    1.05e15         4.39
5 IND    7.15e14         2.99

Then we plot

Code

underestimation_rate <- crop_data[["aggr"]] %>%
filter(ctr == "BRA" |ctr == "USA" | ctr == "IND" | ctr == "ARG" | ctr == "CHN") %>%
  mutate(error_perc = (divtrend_obs - 1)*100 - (divtrend_sim - 1)*100) %>%
  mutate(error_perc = ifelse(is.infinite(error_perc), -1000, error_perc)) %>%  # Replace Inf with -1000
  group_by(model, climate_extreme) %>%
  summarise(under_rate = mean(error_perc < 0) * 100)  # Calculate % of values below 0

dat_long <- crop_data[["aggr"]] %>% 
filter(ctr == "BRA" |ctr == "USA" | ctr == "IND" | ctr == "ARG" | ctr == "CHN") %>%
mutate(divtrend_sim_perc = (divtrend_sim - 1)*100, # Transform to percentages
       divtrend_obs_perc = (divtrend_obs - 1)*100
      ) %>%
  dplyr::select(lat, lon, year, divtrend_sim_perc, divtrend_obs_perc, model, climate_extreme) %>% 
  pivot_longer(
    cols = c(divtrend_obs_perc, divtrend_sim_perc),  # Adjust to the actual column names you want to pivot
    names_to = "type",
    values_to = "divtrend"
  ) %>% 
  mutate(
    model = case_when(
      type == "divtrend_obs_perc" ~ "benchmark",  
      type == "divtrend_sim_perc" ~ model
    )
  ) %>% 
  distinct() 


# Calculate the median yield anomaly for the benchmark model for each climate extreme
benchmark_medians <- dat_long %>%
  filter(model == benchmark_name) %>%
  group_by(climate_extreme) %>%
  summarise(median_yield = median(divtrend, na.rm = TRUE))


# Create the raincloud plot
raincloud_perc <- dat_long %>%
    mutate(model = factor(model, levels = order)) %>%  # Order the models using the defined order
   filter(is.finite(divtrend), abs(divtrend) < 1000) %>%  # Filter out Inf values and keep within range
  ggplot(aes(x = divtrend, y = model, fill = model == benchmark_name)) +
  
  
  # Add half-eye plot
  ggdist::stat_halfeye(aes(fill = model == benchmark_name), adjust = 0.7, 
                       width = 2,               # Increase distribution width
                       justification = -0.1,    # Move half-eye plot further from boxplot
                       .width = c(0.3, 0), 
                       point_colour = "black", 
                       alpha = 0.7) +
  
    
  scale_fill_manual(values = c("TRUE" = "#00A087FF",  
                             "FALSE" = "#7E6148FF"))+ 
  
  # Add boxplots (smaller width to give more space to distribution)
  geom_boxplot(width = 0.12, outlier.shape = NA, alpha = 0.7, 
               fill = "white", color = "black", size = 1) +
  
   # Add red vertical line for the benchmark median
  geom_vline(data = benchmark_medians, aes(xintercept = median_yield), 
              color = "black", size = 2) +

  # Add text annotations for the underestimation rate (place them to the right of the plot)
  geom_text(data = underestimation_rate, 
           aes(x = 40, y = model, label = paste0(round(under_rate, 1), "%")),  # Round percentage to 1 decimal
            vjust = -2, size = 5, aspect = 0.2, color = "black") +  # Adjust text size and positioning
 
  # Labels and annotations

  labs(title = "", 
       subtitle = "Top 5 producers: CHN (23.5%), USA (23.1%), BRA (7.14%), ARG(4.39%), IND(2.99%)",
       x = "Yield anomalies in %", 
       y = "") +
  
  # Add facets for different climate extremes
  facet_wrap(~ climate_extreme, scales = "free_x", ncol = 2) +

  # Custom x-axis breaks
  scale_x_continuous(breaks = seq(-100, 100, by = 20)) +
  
   # Truncate the x-axis to focus on the range
  coord_cartesian(xlim = c(-100, 100)) +
  
   theme(legend.position = "none", 
        axis.title.x = element_text(size = 18),          # Increase x-axis title size
        axis.text.x = element_text(size = 18),           # Increase x-axis text size
        axis.text.y = element_text(size = 18),           # Increase y-axis text size
        strip.text = element_text(size = 18))            # Increase facet label size

raincloud_perc

Code

# Save the figure 
ggsave(file.path(base_path, "suppfig14_RaincloudExtremes.png"), raincloud_perc, width = 18, height = 20, dpi = 300)

RICE

First we calculate who the 5 largest producing countries are and also calculate how much tons they produce and what the global share is (percentage)

Code

crop_datageneral[["ri1"]] %>% 
  mutate(prod = yield * total_area) %>%
  group_by(ctr) %>%
  summarise(ctr_prod = sum(prod, na.rm = TRUE), .groups = "drop") %>%
  mutate(global_prod = sum(ctr_prod),
         global_share = (ctr_prod / global_prod) *100) %>%
  slice_max(order_by = ctr_prod, n = 5) %>%
  dplyr::select(ctr, ctr_prod, global_share)

# A tibble: 5 × 3
  ctr   ctr_prod global_share
  <chr>    <dbl>        <dbl>
1 CHN    8.90e14        30.7 
2 IND    4.11e14        14.2 
3 IDN    2.02e14         6.98
4 THA    1.95e14         6.71
5 BRA    1.84e14         6.35

Then we plot

Code

underestimation_rate <- crop_data[["ri1"]] %>%
filter(ctr == "BRA" |ctr == "THA" | ctr == "IND" | ctr == "IDN" | ctr == "CHN") %>%
  mutate(error_perc = (divtrend_obs - 1)*100 - (divtrend_sim - 1)*100) %>%
  mutate(error_perc = ifelse(is.infinite(error_perc), -1000, error_perc)) %>%  # Replace Inf with -1000
  group_by(model, climate_extreme) %>%
  summarise(under_rate = mean(error_perc < 0) * 100)  # Calculate % of values below 0

dat_long <- crop_data[["ri1"]] %>% 
filter(ctr == "BRA" |ctr == "THA" | ctr == "IND" | ctr == "IDN" | ctr == "CHN") %>%
mutate(divtrend_sim_perc = (divtrend_sim - 1)*100, # Transform to percentages
       divtrend_obs_perc = (divtrend_obs - 1)*100
      ) %>%
  dplyr::select(lat, lon, year, divtrend_sim_perc, divtrend_obs_perc, model, climate_extreme) %>% 
  pivot_longer(
    cols = c(divtrend_obs_perc, divtrend_sim_perc),  # Adjust to the actual column names you want to pivot
    names_to = "type",
    values_to = "divtrend"
  ) %>% 
  mutate(
    model = case_when(
      type == "divtrend_obs_perc" ~ "benchmark",  
      type == "divtrend_sim_perc" ~ model
    )
  ) %>% 
  distinct() 


# Calculate the median yield anomaly for the benchmark model for each climate extreme
benchmark_medians <- dat_long %>%
  filter(model == benchmark_name) %>%
  group_by(climate_extreme) %>%
  summarise(median_yield = median(divtrend, na.rm = TRUE))


# Create the raincloud plot
raincloud_perc <- dat_long %>%
    mutate(model = factor(model, levels = order)) %>%  # Order the models using the defined order
   filter(is.finite(divtrend), abs(divtrend) < 1000) %>%  # Filter out Inf values and keep within range
  ggplot(aes(x = divtrend, y = model, fill = model == benchmark_name)) +
  
  
  # Add half-eye plot
  ggdist::stat_halfeye(aes(fill = model == benchmark_name), adjust = 0.7, 
                       width = 2,               # Increase distribution width
                       justification = -0.1,    # Move half-eye plot further from boxplot
                       .width = c(0.3, 0), 
                       point_colour = "black", 
                       alpha = 0.7) +
  
    
  scale_fill_manual(values = c("TRUE" = "#00A087FF",  
                             "FALSE" = "#7E6148FF"))+ 
  
  # Add boxplots (smaller width to give more space to distribution)
  geom_boxplot(width = 0.12, outlier.shape = NA, alpha = 0.7, 
               fill = "white", color = "black", size = 1) +
  
   # Add red vertical line for the benchmark median
  geom_vline(data = benchmark_medians, aes(xintercept = median_yield), 
              color = "black", size = 2) +

  # Add text annotations for the underestimation rate (place them to the right of the plot)
  geom_text(data = underestimation_rate, 
           aes(x = 40, y = model, label = paste0(round(under_rate, 1), "%")),  # Round percentage to 1 decimal
            vjust = -2, size = 5, aspect = 0.2, color = "black") +  # Adjust text size and positioning
 
  # Labels and annotations

  labs(title = "", 
       subtitle = "Top 5 producers: CHN (30.7%), IND (14.2%), IDN (6.98%), THA(6.71%), BRA(6.35%)",
       x = "Yield anomalies in %", 
       y = "") +
  
  # Add facets for different climate extremes
  facet_wrap(~ climate_extreme, scales = "free_x", ncol = 2) +

  # Custom x-axis breaks
  scale_x_continuous(breaks = seq(-100, 100, by = 20)) +
  
   # Truncate the x-axis to focus on the range
  coord_cartesian(xlim = c(-100, 100)) +
  
   theme(legend.position = "none", 
        axis.title.x = element_text(size = 18),          # Increase x-axis title size
        axis.text.x = element_text(size = 18),           # Increase x-axis text size
        axis.text.y = element_text(size = 18),           # Increase y-axis text size
        strip.text = element_text(size = 18))            # Increase facet label size

raincloud_perc

Code

# Save the figure 
ggsave(file.path(base_path, "suppfig15_RaincloudExtremes.png"), raincloud_perc, width = 18, height = 20, dpi = 300)

WHEAT

First we calculate who the 5 largest producing countries are and also calculate how much tons they produce and what the global share is (percentage)

Code

crop_datageneral[["wwh"]] %>% 
  mutate(prod = yield * total_area) %>%
  group_by(ctr) %>%
  summarise(ctr_prod = sum(prod, na.rm = TRUE), .groups = "drop") %>%
  mutate(global_prod = sum(ctr_prod),
         global_share = (ctr_prod / global_prod) *100) %>%
  slice_max(order_by = ctr_prod, n = 5) %>%
  dplyr::select(ctr, ctr_prod, global_share)

# A tibble: 5 × 3
  ctr   ctr_prod global_share
  <chr>    <dbl>        <dbl>
1 USA    4.15e14        14.1 
2 FRA    3.75e14        12.7 
3 CHN    3.35e14        11.4 
4 DEU    2.51e14         8.53
5 TUR    1.76e14         5.97

Then we plot

Code

underestimation_rate <- crop_data[["wwh"]] %>%
filter(ctr == "FRA" |ctr == "USA" | ctr == "DEU" | ctr == "TUR" | ctr == "CHN") %>%
  mutate(error_perc = (divtrend_obs - 1)*100 - (divtrend_sim - 1)*100) %>%
  mutate(error_perc = ifelse(is.infinite(error_perc), -1000, error_perc)) %>%  # Replace Inf with -1000
  group_by(model, climate_extreme) %>%
  summarise(under_rate = mean(error_perc < 0) * 100)  # Calculate % of values below 0

dat_long <- crop_data[["wwh"]] %>% 
filter(ctr == "FRA" |ctr == "USA" | ctr == "DEU" | ctr == "TUR" | ctr == "CHN") %>%
mutate(divtrend_sim_perc = (divtrend_sim - 1)*100, # Transform to percentages
       divtrend_obs_perc = (divtrend_obs - 1)*100
      ) %>%
  dplyr::select(lat, lon, year, divtrend_sim_perc, divtrend_obs_perc, model, climate_extreme) %>% 
  pivot_longer(
    cols = c(divtrend_obs_perc, divtrend_sim_perc),  # Adjust to the actual column names you want to pivot
    names_to = "type",
    values_to = "divtrend"
  ) %>% 
  mutate(
    model = case_when(
      type == "divtrend_obs_perc" ~ "benchmark",  
      type == "divtrend_sim_perc" ~ model
    )
  ) %>% 
  distinct() 


# Calculate the median yield anomaly for the benchmark model for each climate extreme
benchmark_medians <- dat_long %>%
  filter(model == benchmark_name) %>%
  group_by(climate_extreme) %>%
  summarise(median_yield = median(divtrend, na.rm = TRUE))


# Create the raincloud plot
raincloud_perc <- dat_long %>%
    mutate(model = factor(model, levels = order)) %>%  # Order the models using the defined order
   filter(is.finite(divtrend), abs(divtrend) < 1000) %>%  # Filter out Inf values and keep within range
  ggplot(aes(x = divtrend, y = model, fill = model == benchmark_name)) +
  
  
  # Add half-eye plot
  ggdist::stat_halfeye(aes(fill = model == benchmark_name), adjust = 0.7, 
                       width = 2,               # Increase distribution width
                       justification = -0.1,    # Move half-eye plot further from boxplot
                       .width = c(0.3, 0), 
                       point_colour = "black", 
                       alpha = 0.7) +
  
    
  scale_fill_manual(values = c("TRUE" = "#00A087FF",  
                             "FALSE" = "#7E6148FF"))+ 
  
  # Add boxplots (smaller width to give more space to distribution)
  geom_boxplot(width = 0.12, outlier.shape = NA, alpha = 0.7, 
               fill = "white", color = "black", size = 1) +
  
   # Add red vertical line for the benchmark median
  geom_vline(data = benchmark_medians, aes(xintercept = median_yield), 
              color = "black", size = 2) +

  # Add text annotations for the underestimation rate (place them to the right of the plot)
  geom_text(data = underestimation_rate, 
           aes(x = 40, y = model, label = paste0(round(under_rate, 1), "%")),  # Round percentage to 1 decimal
            vjust = -2, size = 5, aspect = 0.2, color = "black") +  # Adjust text size and positioning
 
  # Labels and annotations

  labs(title = "", 
       subtitle = "Top 5 producers: USA (14.1%), FRA (12.7%), CHN (11.4%), DEU(8.53%), TUR(5.97%)",
       x = "Yield anomalies in %", 
       y = "") +
  
  # Add facets for different climate extremes
  facet_wrap(~ climate_extreme, scales = "free_x", ncol = 2) +

  # Custom x-axis breaks
  scale_x_continuous(breaks = seq(-100, 100, by = 20)) +
  
   # Truncate the x-axis to focus on the range
  coord_cartesian(xlim = c(-100, 100)) +
  
   theme(legend.position = "none", 
        axis.title.x = element_text(size = 18),          # Increase x-axis title size
        axis.text.x = element_text(size = 18),           # Increase x-axis text size
        axis.text.y = element_text(size = 18),           # Increase y-axis text size
        strip.text = element_text(size = 18))            # Increase facet label size

raincloud_perc

Code

# Save the figure 
ggsave(file.path(base_path, "suppfig16_RaincloudExtremes.png"), raincloud_perc, width = 18, height = 20, dpi = 300)

SOY

First we calculate who the 5 largest producing countries are and also calculate how much tons they produce and what the global share is (percentage)

Code

crop_datageneral[["soy"]] %>% 
  mutate(prod = yield * total_area) %>%
  group_by(ctr) %>%
  summarise(ctr_prod = sum(prod, na.rm = TRUE), .groups = "drop") %>%
  mutate(global_prod = sum(ctr_prod),
         global_share = (ctr_prod / global_prod) *100) %>%
  slice_max(order_by = ctr_prod, n = 5) %>%
  dplyr::select(ctr, ctr_prod, global_share)

# A tibble: 5 × 3
  ctr   ctr_prod global_share
  <chr>    <dbl>        <dbl>
1 USA    5.92e14        41.4 
2 BRA    4.04e14        28.2 
3 ARG    2.40e14        16.8 
4 CHN    1.11e14         7.78
5 IND    1.96e13         1.37

Then we plot

Code

underestimation_rate <- crop_data[["soy"]] %>%
filter(ctr == "BRA" |ctr == "USA" | ctr == "IND" | ctr == "ARG" | ctr == "CHN") %>%
  mutate(error_perc = (divtrend_obs - 1)*100 - (divtrend_sim - 1)*100) %>%
  mutate(error_perc = ifelse(is.infinite(error_perc), -1000, error_perc)) %>%  # Replace Inf with -1000
  group_by(model, climate_extreme) %>%
  summarise(under_rate = mean(error_perc < 0) * 100)  # Calculate % of values below 0

dat_long <- crop_data[["soy"]] %>% 
filter(ctr == "BRA" |ctr == "USA" | ctr == "IND" | ctr == "ARG" | ctr == "CHN") %>%
mutate(divtrend_sim_perc = (divtrend_sim - 1)*100, # Transform to percentages
       divtrend_obs_perc = (divtrend_obs - 1)*100
      ) %>%
  dplyr::select(lat, lon, year, divtrend_sim_perc, divtrend_obs_perc, model, climate_extreme) %>% 
  pivot_longer(
    cols = c(divtrend_obs_perc, divtrend_sim_perc),  # Adjust to the actual column names you want to pivot
    names_to = "type",
    values_to = "divtrend"
  ) %>% 
  mutate(
    model = case_when(
      type == "divtrend_obs_perc" ~ "benchmark",  
      type == "divtrend_sim_perc" ~ model
    )
  ) %>% 
  distinct() 


# Calculate the median yield anomaly for the benchmark model for each climate extreme
benchmark_medians <- dat_long %>%
  filter(model == benchmark_name) %>%
  group_by(climate_extreme) %>%
  summarise(median_yield = median(divtrend, na.rm = TRUE))


# Create the raincloud plot
raincloud_perc <- dat_long %>%
    mutate(model = factor(model, levels = order)) %>%  # Order the models using the defined order
   filter(is.finite(divtrend), abs(divtrend) < 1000) %>%  # Filter out Inf values and keep within range
  ggplot(aes(x = divtrend, y = model, fill = model == benchmark_name)) +
  
  
  # Add half-eye plot
  ggdist::stat_halfeye(aes(fill = model == benchmark_name), adjust = 0.7, 
                       width = 2,               # Increase distribution width
                       justification = -0.1,    # Move half-eye plot further from boxplot
                       .width = c(0.3, 0), 
                       point_colour = "black", 
                       alpha = 0.7) +
  
    
  scale_fill_manual(values = c("TRUE" = "#00A087FF",  
                             "FALSE" = "#7E6148FF"))+ 
  
  # Add boxplots (smaller width to give more space to distribution)
  geom_boxplot(width = 0.12, outlier.shape = NA, alpha = 0.7, 
               fill = "white", color = "black", size = 1) +
  
   # Add red vertical line for the benchmark median
  geom_vline(data = benchmark_medians, aes(xintercept = median_yield), 
              color = "black", size = 2) +

  # Add text annotations for the underestimation rate (place them to the right of the plot)
  geom_text(data = underestimation_rate, 
           aes(x = 40, y = model, label = paste0(round(under_rate, 1), "%")),  # Round percentage to 1 decimal
            vjust = -2, size = 5, aspect = 0.2, color = "black") +  # Adjust text size and positioning
 
  # Labels and annotations

  labs(title = "", 
       subtitle = "Top 5 producers: USA (41.4%), BRA (28.2%), ARG (16.8%), CHN(7.78%), IND(1.37%)",
       x = "Yield anomalies in %", 
       y = "") +
  
  # Add facets for different climate extremes
  facet_wrap(~ climate_extreme, scales = "free_x", ncol = 2) +

  # Custom x-axis breaks
  scale_x_continuous(breaks = seq(-100, 100, by = 20)) +
  
   # Truncate the x-axis to focus on the range
  coord_cartesian(xlim = c(-100, 100)) +
  
   theme(legend.position = "none", 
        axis.title.x = element_text(size = 18),          # Increase x-axis title size
        axis.text.x = element_text(size = 18),           # Increase x-axis text size
        axis.text.y = element_text(size = 18),           # Increase y-axis text size
        strip.text = element_text(size = 18))            # Increase facet label size

raincloud_perc

Code

# Save the figure 
ggsave(file.path(base_path, "suppfig17_RaincloudExtremes.png"), raincloud_perc, width = 18, height = 20, dpi = 300)

MAIZE

First we calculate who the 5 largest producing countries are and also calculate how much tons they produce and what the global share is (percentage)

Code

crop_datageneral[["mai"]] %>% 
  mutate(prod = yield * total_area) %>%
  group_by(ctr) %>%
  summarise(ctr_prod = sum(prod, na.rm = TRUE), .groups = "drop") %>%
  mutate(global_prod = sum(ctr_prod),
         global_share = (ctr_prod / global_prod) *100) %>%
  slice_max(order_by = ctr_prod, n = 5) %>%
  dplyr::select(ctr, ctr_prod, global_share)

# A tibble: 5 × 3
  ctr   ctr_prod global_share
  <chr>    <dbl>        <dbl>
1 USA    2.30e15        40.6 
2 CHN    9.30e14        16.4 
3 BRA    4.49e14         7.91
4 ARG    2.47e14         4.35
5 FRA    2.41e14         4.24

Then we plot

Code

underestimation_rate <- crop_data[["mai"]] %>%
filter(ctr == "BRA" |ctr == "USA" | ctr == "FRA" | ctr == "ARG" | ctr == "CHN") %>%
  mutate(error_perc = (divtrend_obs - 1)*100 - (divtrend_sim - 1)*100) %>%
  mutate(error_perc = ifelse(is.infinite(error_perc), -1000, error_perc)) %>%  # Replace Inf with -1000
  group_by(model, climate_extreme) %>%
  summarise(under_rate = mean(error_perc < 0) * 100)  # Calculate % of values below 0

dat_long <- crop_data[["mai"]] %>% 
filter(ctr == "BRA" |ctr == "USA" | ctr == "FRA" | ctr == "ARG" | ctr == "CHN") %>%
mutate(divtrend_sim_perc = (divtrend_sim - 1)*100, # Transform to percentages
       divtrend_obs_perc = (divtrend_obs - 1)*100
      ) %>%
  dplyr::select(lat, lon, year, divtrend_sim_perc, divtrend_obs_perc, model, climate_extreme) %>% 
  pivot_longer(
    cols = c(divtrend_obs_perc, divtrend_sim_perc),  # Adjust to the actual column names you want to pivot
    names_to = "type",
    values_to = "divtrend"
  ) %>% 
  mutate(
    model = case_when(
      type == "divtrend_obs_perc" ~ "benchmark",  
      type == "divtrend_sim_perc" ~ model
    )
  ) %>% 
  distinct() 


# Calculate the median yield anomaly for the benchmark model for each climate extreme
benchmark_medians <- dat_long %>%
  filter(model == benchmark_name) %>%
  group_by(climate_extreme) %>%
  summarise(median_yield = median(divtrend, na.rm = TRUE))


# Create the raincloud plot
raincloud_perc <- dat_long %>%
    mutate(model = factor(model, levels = order)) %>%  # Order the models using the defined order
   filter(is.finite(divtrend), abs(divtrend) < 1000) %>%  # Filter out Inf values and keep within range
  ggplot(aes(x = divtrend, y = model, fill = model == benchmark_name)) +
  
  
  # Add half-eye plot
  ggdist::stat_halfeye(aes(fill = model == benchmark_name), adjust = 0.7, 
                       width = 2,               # Increase distribution width
                       justification = -0.1,    # Move half-eye plot further from boxplot
                       .width = c(0.3, 0), 
                       point_colour = "black", 
                       alpha = 0.7) +
  
    
  scale_fill_manual(values = c("TRUE" = "#00A087FF",  
                             "FALSE" = "#7E6148FF"))+ 
  
  # Add boxplots (smaller width to give more space to distribution)
  geom_boxplot(width = 0.12, outlier.shape = NA, alpha = 0.7, 
               fill = "white", color = "black", size = 1) +
  
   # Add red vertical line for the benchmark median
  geom_vline(data = benchmark_medians, aes(xintercept = median_yield), 
              color = "black", size = 2) +

  # Add text annotations for the underestimation rate (place them to the right of the plot)
  geom_text(data = underestimation_rate, 
           aes(x = 40, y = model, label = paste0(round(under_rate, 1), "%")),  # Round percentage to 1 decimal
            vjust = -2, size = 5, aspect = 0.2, color = "black") +  # Adjust text size and positioning
 
  # Labels and annotations

  labs(title = "", 
       subtitle = "Top 5 producers: USA (40.6%), CHN (16.4%), BRA (7.91%), ARG(4.35%), FRA(4.24%)",
       x = "Yield anomalies in %", 
       y = "") +
  
  # Add facets for different climate extremes
  facet_wrap(~ climate_extreme, scales = "free_x", ncol = 2) +

  # Custom x-axis breaks
  scale_x_continuous(breaks = seq(-100, 100, by = 20)) +
  
   # Truncate the x-axis to focus on the range
  coord_cartesian(xlim = c(-100, 100)) +
  
   theme(legend.position = "none", 
        axis.title.x = element_text(size = 18),          # Increase x-axis title size
        axis.text.x = element_text(size = 18),           # Increase x-axis text size
        axis.text.y = element_text(size = 18),           # Increase y-axis text size
        strip.text = element_text(size = 18))            # Increase facet label size

raincloud_perc

Code

# Save the figure 
ggsave(file.path(base_path, "suppfig18_RaincloudExtremes.png"), raincloud_perc, width = 18, height = 20, dpi = 300)

Supplementary figures 19-23 – rainclouds for absolute detrending

AGGREGATED

Code

dat_long <- crop_data[["aggr"]] %>% 
  dplyr::select(lat, lon, year, difftrend_sim, difftrend_obs, model, climate_extreme) %>%
  pivot_longer(
    cols = c(difftrend_obs, difftrend_sim),  # Adjust to the actual column names you want to pivot
    names_to = "type",
    values_to = "difftrend"
  ) %>% 
  mutate(
    model = case_when(
      type == "difftrend_obs" ~ "benchmark",  
      type == "difftrend_sim" ~ model
    )
  ) %>% 
  distinct() 

underestimation_rate <- crop_data[["aggr"]] %>%
  mutate(error = difftrend_obs - difftrend_sim) %>%
  group_by(model, climate_extreme) %>%
  summarise(under_rate = mean(error < 0) * 100)  # Calculate % of values below 0

benchmark_medians <- dat_long %>%
  filter(model == "benchmark") %>%
  group_by(climate_extreme) %>%
  summarise(median_yield = median(difftrend, na.rm = TRUE))

# Create the raincloud plot
raincloud_perc <- dat_long %>%
    mutate(model = factor(model, levels = order)) %>%  # Order the models using the defined order
  ggplot(aes(x = difftrend, y = model, fill = model == benchmark_name)) +
  
  
  # Add half-eye plot
  ggdist::stat_halfeye(aes(fill = model == benchmark_name), adjust = 0.7, 
                       width = 0.7,               # Increase distribution width
                       justification = -0.1,    # Move half-eye plot further from boxplot
                       .width = c(0.3, 0), 
                       point_colour = "black", 
                       alpha = 0.7) +
  
    
  scale_fill_manual(values = c("TRUE" = "#00A087FF",  
                             "FALSE" = "#7E6148FF"))+ 
  
  # Add boxplots (smaller width to give more space to distribution)
  geom_boxplot(width = 0.12, outlier.shape = NA, alpha = 0.7, 
               fill = "white", color = "black", size = 1) +
  
   # Add black vertical line for the benchmark median
  geom_vline(data = benchmark_medians, aes(xintercept = median_yield), 
              color = "black", size = 2) +
 
  # Add text annotations for the underestimation rate (place them to the right of the plot)
  geom_text(data = underestimation_rate, 
           aes(x = 40, y = model, label = paste0(round(under_rate, 1), "%")),  # Round percentage to 1 decimal
            vjust = -2, size = 5, color = "black") +  # Adjust text size and positioning
  
  # Labels and annotations
  labs(title = "", 
       x = "Yield anomalies in tons/ha", 
       y = "") +
  
  # Add facets for different climate extremes
  facet_wrap(~ climate_extreme, scales = "free_x", ncol = 2) +
  
 # Custom x-axis breaks
  scale_x_continuous(breaks = seq(-3, 3, by = 1)) +
  
   # Truncate the x-axis to focus on the range
  coord_cartesian(xlim = c(-3, 3)) +
  
   theme(legend.position = "none", 
        axis.title.x = element_text(size = 18),          # Increase x-axis title size
        axis.text.x = element_text(size = 18),           # Increase x-axis text size
        axis.text.y = element_text(size = 18),           # Increase y-axis text size
        strip.text = element_text(size = 18))            # Increase facet label size

raincloud_perc

Code

# Save the figure 
ggsave(file.path(base_path, "suppfig19_RaincloudExtremes.png"), raincloud_perc, width = 18, height = 20, dpi = 300)

RICE

Code

dat_long <- crop_data[["ri1"]] %>% 
  dplyr::select(lat, lon, year, difftrend_sim, difftrend_obs, model, climate_extreme) %>%
  pivot_longer(
    cols = c(difftrend_obs, difftrend_sim),  # Adjust to the actual column names you want to pivot
    names_to = "type",
    values_to = "difftrend"
  ) %>% 
  mutate(
    model = case_when(
      type == "difftrend_obs" ~ "benchmark",  
      type == "difftrend_sim" ~ model
    )
  ) %>% 
  distinct() 

underestimation_rate <- crop_data[["ri1"]] %>%
  mutate(error = difftrend_obs - difftrend_sim) %>%
  group_by(model, climate_extreme) %>%
  summarise(under_rate = mean(error < 0) * 100)  # Calculate % of values below 0

benchmark_medians <- dat_long %>%
  filter(model == "benchmark") %>%
  group_by(climate_extreme) %>%
  summarise(median_yield = median(difftrend, na.rm = TRUE))

# Create the raincloud plot
raincloud_perc <- dat_long %>%
    mutate(model = factor(model, levels = order)) %>%  # Order the models using the defined order
  ggplot(aes(x = difftrend, y = model, fill = model == benchmark_name)) +
  
  
  # Add half-eye plot
  ggdist::stat_halfeye(aes(fill = model == benchmark_name), adjust = 0.7, 
                       width = 0.7,               # Increase distribution width
                       justification = -0.1,    # Move half-eye plot further from boxplot
                       .width = c(0.3, 0), 
                       point_colour = "black", 
                       alpha = 0.7) +
  
    
  scale_fill_manual(values = c("TRUE" = "#00A087FF",  
                             "FALSE" = "#7E6148FF"))+ 
  
  # Add boxplots (smaller width to give more space to distribution)
  geom_boxplot(width = 0.12, outlier.shape = NA, alpha = 0.7, 
               fill = "white", color = "black", size = 1) +
  
   # Add black vertical line for the benchmark median
  geom_vline(data = benchmark_medians, aes(xintercept = median_yield), 
              color = "black", size = 2) +
 
  # Add text annotations for the underestimation rate (place them to the right of the plot)
  geom_text(data = underestimation_rate, 
           aes(x = 40, y = model, label = paste0(round(under_rate, 1), "%")),  # Round percentage to 1 decimal
            vjust = -2, size = 5, color = "black") +  # Adjust text size and positioning
  
  # Labels and annotations
  labs(title = "", 
       x = "Yield anomalies in tons/ha", 
       y = "") +
  
  # Add facets for different climate extremes
  facet_wrap(~ climate_extreme, scales = "free_x", ncol = 2) +
  
 # Custom x-axis breaks
  scale_x_continuous(breaks = seq(-3, 3, by = 1)) +
  
   # Truncate the x-axis to focus on the range
  coord_cartesian(xlim = c(-3, 3)) +
  
   theme(legend.position = "none", 
        axis.title.x = element_text(size = 18),          # Increase x-axis title size
        axis.text.x = element_text(size = 18),           # Increase x-axis text size
        axis.text.y = element_text(size = 18),           # Increase y-axis text size
        strip.text = element_text(size = 18))            # Increase facet label size

raincloud_perc

Code

# Save the figure 
ggsave(file.path(base_path, "suppfig20_RaincloudExtremes.png"), raincloud_perc, width = 18, height = 20, dpi = 300)

WHEAT

Code

dat_long <- crop_data[["wwh"]] %>% 
  dplyr::select(lat, lon, year, difftrend_sim, difftrend_obs, model, climate_extreme) %>%
  pivot_longer(
    cols = c(difftrend_obs, difftrend_sim),  # Adjust to the actual column names you want to pivot
    names_to = "type",
    values_to = "difftrend"
  ) %>% 
  mutate(
    model = case_when(
      type == "difftrend_obs" ~ "benchmark",  
      type == "difftrend_sim" ~ model
    )
  ) %>% 
  distinct() 

underestimation_rate <- crop_data[["wwh"]] %>%
  mutate(error = difftrend_obs - difftrend_sim) %>%
  group_by(model, climate_extreme) %>%
  summarise(under_rate = mean(error < 0) * 100)  # Calculate % of values below 0

benchmark_medians <- dat_long %>%
  filter(model == "benchmark") %>%
  group_by(climate_extreme) %>%
  summarise(median_yield = median(difftrend, na.rm = TRUE))

# Create the raincloud plot
raincloud_perc <- dat_long %>%
    mutate(model = factor(model, levels = order)) %>%  # Order the models using the defined order
  ggplot(aes(x = difftrend, y = model, fill = model == benchmark_name)) +
  
  
  # Add half-eye plot
  ggdist::stat_halfeye(aes(fill = model == benchmark_name), adjust = 0.7, 
                       width = 0.7,               # Increase distribution width
                       justification = -0.1,    # Move half-eye plot further from boxplot
                       .width = c(0.3, 0), 
                       point_colour = "black", 
                       alpha = 0.7) +
  
    
  scale_fill_manual(values = c("TRUE" = "#00A087FF",  
                             "FALSE" = "#7E6148FF"))+ 
  
  # Add boxplots (smaller width to give more space to distribution)
  geom_boxplot(width = 0.12, outlier.shape = NA, alpha = 0.7, 
               fill = "white", color = "black", size = 1) +
  
   # Add black vertical line for the benchmark median
  geom_vline(data = benchmark_medians, aes(xintercept = median_yield), 
              color = "black", size = 2) +
 
  # Add text annotations for the underestimation rate (place them to the right of the plot)
  geom_text(data = underestimation_rate, 
           aes(x = 40, y = model, label = paste0(round(under_rate, 1), "%")),  # Round percentage to 1 decimal
            vjust = -2, size = 5, color = "black") +  # Adjust text size and positioning
  
  # Labels and annotations
  labs(title = "", 
       x = "Yield anomalies in tons/ha", 
       y = "") +
  
  # Add facets for different climate extremes
  facet_wrap(~ climate_extreme, scales = "free_x", ncol = 2) +
  
 # Custom x-axis breaks
  scale_x_continuous(breaks = seq(-3, 3, by = 1)) +
  
   # Truncate the x-axis to focus on the range
  coord_cartesian(xlim = c(-3, 3)) +
  
   theme(legend.position = "none", 
        axis.title.x = element_text(size = 18),          # Increase x-axis title size
        axis.text.x = element_text(size = 18),           # Increase x-axis text size
        axis.text.y = element_text(size = 18),           # Increase y-axis text size
        strip.text = element_text(size = 18))            # Increase facet label size

raincloud_perc

Code

# Save the figure 
ggsave(file.path(base_path, "suppfig21_RaincloudExtremes.png"), raincloud_perc, width = 18, height = 20, dpi = 300)

SOY

Code

dat_long <- crop_data[["soy"]] %>% 
  dplyr::select(lat, lon, year, difftrend_sim, difftrend_obs, model, climate_extreme) %>%
  pivot_longer(
    cols = c(difftrend_obs, difftrend_sim),  # Adjust to the actual column names you want to pivot
    names_to = "type",
    values_to = "difftrend"
  ) %>% 
  mutate(
    model = case_when(
      type == "difftrend_obs" ~ "benchmark",  
      type == "difftrend_sim" ~ model
    )
  ) %>% 
  distinct() 

underestimation_rate <- crop_data[["soy"]] %>%
  mutate(error = difftrend_obs - difftrend_sim) %>%
  group_by(model, climate_extreme) %>%
  summarise(under_rate = mean(error < 0) * 100)  # Calculate % of values below 0

benchmark_medians <- dat_long %>%
  filter(model == "benchmark") %>%
  group_by(climate_extreme) %>%
  summarise(median_yield = median(difftrend, na.rm = TRUE))

# Create the raincloud plot
raincloud_perc <- dat_long %>%
    mutate(model = factor(model, levels = order)) %>%  # Order the models using the defined order
  ggplot(aes(x = difftrend, y = model, fill = model == benchmark_name)) +
  
  
  # Add half-eye plot
  ggdist::stat_halfeye(aes(fill = model == benchmark_name), adjust = 0.7, 
                       width = 0.7,               # Increase distribution width
                       justification = -0.1,    # Move half-eye plot further from boxplot
                       .width = c(0.3, 0), 
                       point_colour = "black", 
                       alpha = 0.7) +
  
    
  scale_fill_manual(values = c("TRUE" = "#00A087FF",  
                             "FALSE" = "#7E6148FF"))+ 
  
  # Add boxplots (smaller width to give more space to distribution)
  geom_boxplot(width = 0.12, outlier.shape = NA, alpha = 0.7, 
               fill = "white", color = "black", size = 1) +
  
   # Add black vertical line for the benchmark median
  geom_vline(data = benchmark_medians, aes(xintercept = median_yield), 
              color = "black", size = 2) +
 
  # Add text annotations for the underestimation rate (place them to the right of the plot)
  geom_text(data = underestimation_rate, 
           aes(x = 40, y = model, label = paste0(round(under_rate, 1), "%")),  # Round percentage to 1 decimal
            vjust = -2, size = 5, color = "black") +  # Adjust text size and positioning
  
  # Labels and annotations
  labs(title = "", 
       x = "Yield anomalies in tons/ha", 
       y = "") +
  
  # Add facets for different climate extremes
  facet_wrap(~ climate_extreme, scales = "free_x", ncol = 2) +
  
 # Custom x-axis breaks
  scale_x_continuous(breaks = seq(-3, 3, by = 1)) +
  
   # Truncate the x-axis to focus on the range
  coord_cartesian(xlim = c(-3, 3)) +
  
   theme(legend.position = "none", 
        axis.title.x = element_text(size = 18),          # Increase x-axis title size
        axis.text.x = element_text(size = 18),           # Increase x-axis text size
        axis.text.y = element_text(size = 18),           # Increase y-axis text size
        strip.text = element_text(size = 18))            # Increase facet label size

raincloud_perc

Code

# Save the figure 
ggsave(file.path(base_path, "suppfig22_RaincloudExtremes.png"), raincloud_perc, width = 18, height = 20, dpi = 300)

MAIZE

Code

dat_long <- crop_data[["mai"]] %>% 
  dplyr::select(lat, lon, year, difftrend_sim, difftrend_obs, model, climate_extreme) %>%
  pivot_longer(
    cols = c(difftrend_obs, difftrend_sim),  # Adjust to the actual column names you want to pivot
    names_to = "type",
    values_to = "difftrend"
  ) %>% 
  mutate(
    model = case_when(
      type == "difftrend_obs" ~ "benchmark",  
      type == "difftrend_sim" ~ model
    )
  ) %>% 
  distinct() 

underestimation_rate <- crop_data[["mai"]] %>%
  mutate(error = difftrend_obs - difftrend_sim) %>%
  group_by(model, climate_extreme) %>%
  summarise(under_rate = mean(error < 0) * 100)  # Calculate % of values below 0

benchmark_medians <- dat_long %>%
  filter(model == "benchmark") %>%
  group_by(climate_extreme) %>%
  summarise(median_yield = median(difftrend, na.rm = TRUE))

# Create the raincloud plot
raincloud_perc <- dat_long %>%
    mutate(model = factor(model, levels = order)) %>%  # Order the models using the defined order
  ggplot(aes(x = difftrend, y = model, fill = model == benchmark_name)) +
  
  
  # Add half-eye plot
  ggdist::stat_halfeye(aes(fill = model == benchmark_name), adjust = 0.7, 
                       width = 0.7,               # Increase distribution width
                       justification = -0.1,    # Move half-eye plot further from boxplot
                       .width = c(0.3, 0), 
                       point_colour = "black", 
                       alpha = 0.7) +
  
    
  scale_fill_manual(values = c("TRUE" = "#00A087FF",  
                             "FALSE" = "#7E6148FF"))+ 
  
  # Add boxplots (smaller width to give more space to distribution)
  geom_boxplot(width = 0.12, outlier.shape = NA, alpha = 0.7, 
               fill = "white", color = "black", size = 1) +
  
   # Add black vertical line for the benchmark median
  geom_vline(data = benchmark_medians, aes(xintercept = median_yield), 
              color = "black", size = 2) +
 
  # Add text annotations for the underestimation rate (place them to the right of the plot)
  geom_text(data = underestimation_rate, 
           aes(x = 40, y = model, label = paste0(round(under_rate, 1), "%")),  # Round percentage to 1 decimal
            vjust = -2, size = 5, color = "black") +  # Adjust text size and positioning
  
  # Labels and annotations
  labs(title = "", 
       x = "Yield anomalies in tons/ha", 
       y = "") +
  
  # Add facets for different climate extremes
  facet_wrap(~ climate_extreme, scales = "free_x", ncol = 2) +
  
 # Custom x-axis breaks
  scale_x_continuous(breaks = seq(-3, 3, by = 1)) +
  
   # Truncate the x-axis to focus on the range
  coord_cartesian(xlim = c(-3, 3)) +
  
   theme(legend.position = "none", 
        axis.title.x = element_text(size = 18),          # Increase x-axis title size
        axis.text.x = element_text(size = 18),           # Increase x-axis text size
        axis.text.y = element_text(size = 18),           # Increase y-axis text size
        strip.text = element_text(size = 18))            # Increase facet label size

raincloud_perc

Code

# Save the figure 
ggsave(file.path(base_path, "suppfig23_RaincloudExtremes.png"), raincloud_perc, width = 18, height = 20, dpi = 300)

Supplementary figures 24-28 – rainclouds for relative model error

Instead of separating benchmark from simulation data we can also create the raincloud plots for the model errors (benchmark - simulation data points)

AGGREGATED

Code

underestimation_rate <- crop_data[["aggr"]] %>%
mutate(error_perc = (divtrend_obs - 1)*100 - (divtrend_sim - 1)*100) %>%
  mutate(error_perc = ifelse(is.infinite(error_perc), -1000, error_perc)) %>%  # Replace Inf with -1000
  group_by(model, climate_extreme) %>%
  summarise(under_rate = mean(error_perc < 0) * 100)  # Calculate % of values below 0


# Create the raincloud plot with Inf values filtered out
raincloud_perc <- crop_data[["aggr"]] %>%
mutate(error_perc = (divtrend_obs - 1)*100 - (divtrend_sim - 1)*100) %>%
  filter(is.finite(error_perc), abs(error_perc) < 1000) %>%  # Filter out Inf values and keep within range
  mutate(model = factor(model, levels = order)) %>%  # Order the models using the defined order
  ggplot(aes(x = error_perc, y = model)) +
  coord_trans(x = scales::modulus_trans(p = 0.5)) +  # Adjust p for stronger compression
  
  # Add half-eye plot with fill based on whether the error is negative or positive
  ggdist::stat_halfeye(aes(fill = stat(x < 0)),  # Logical fill based on error sign
                       adjust = 0.7, 
                       width = 0.7,               # Increase distribution width
                       justification = -0.1,    # Move half-eye plot further from boxplot
                       .width = c(0.3, 0), 
                       point_colour = "black", 
                       alpha = 0.7) +
  
  # Add boxplots (smaller width to give more space to distribution)
  geom_boxplot(width = 0.12, outlier.shape = NA, alpha = 0.7, 
               fill = "white", color = "black", size = 1) +
  
  # Add a vertical line at 0 (benchmark median)
  geom_vline(aes(xintercept = 0), color = "black", size = 2) +
  
  # Custom color scale where negative values are red and positive values are blue
  scale_fill_manual(values = c("grey", "#DC0000FF") )+
                      
    # Add text annotations for the underestimation rate (place them to the right of the plot)
  geom_text(data = underestimation_rate, 
           aes(x = -50.25, y = model, label = paste0(round(under_rate, 1), "%")),  # Round percentage to 1 decimal
            vjust = -2, size = 5, aspect = 0.2, color = "#DC0000FF") +  # Adjust text size and positioning
  
  # Labels and annotation
  labs(x = "Model Error (benchmark - simulated) in % yield anomalies", y = "") +
  
  # Add facets for different climate extremes
  facet_wrap(~ climate_extreme, scales = "free_x", ncol = 2) +
  
  # Custom x-axis breaks
  scale_x_continuous(breaks = seq(-100, 100, by = 20)) +
  
  # Truncate the x-axis to focus on the range
  coord_cartesian(xlim = c(-100, 100)) +

  theme(legend.position = "none", 
        axis.title.x = element_text(size = 18),          # Increase x-axis title size
        axis.text.x = element_text(size = 18),           # Increase x-axis text size
        axis.text.y = element_text(size = 18),           # Increase y-axis text size
        strip.text = element_text(size = 18))            # Increase facet label size

raincloud_perc

Code

ggsave(file.path(base_path, "suppfig24_RaincloudExtremes.png"), raincloud_perc, width = 18, height = 20, dpi = 300)

RICE

Code

underestimation_rate <- crop_data[["ri1"]] %>%
mutate(error_perc = (divtrend_obs - 1)*100 - (divtrend_sim - 1)*100) %>%
  mutate(error_perc = ifelse(is.infinite(error_perc), -1000, error_perc)) %>%  # Replace Inf with -1000
  group_by(model, climate_extreme) %>%
  summarise(under_rate = mean(error_perc < 0) * 100)  # Calculate % of values below 0


# Create the raincloud plot with Inf values filtered out
raincloud_perc <- crop_data[["ri1"]] %>%
mutate(error_perc = (divtrend_obs - 1)*100 - (divtrend_sim - 1)*100) %>%
  filter(is.finite(error_perc), abs(error_perc) < 1000) %>%  # Filter out Inf values and keep within range
  mutate(model = factor(model, levels = order)) %>%  # Order the models using the defined order
  ggplot(aes(x = error_perc, y = model)) +
  coord_trans(x = scales::modulus_trans(p = 0.5)) +  # Adjust p for stronger compression
  
  # Add half-eye plot with fill based on whether the error is negative or positive
  ggdist::stat_halfeye(aes(fill = stat(x < 0)),  # Logical fill based on error sign
                       adjust = 0.7, 
                       width = 0.7,               # Increase distribution width
                       justification = -0.1,    # Move half-eye plot further from boxplot
                       .width = c(0.3, 0), 
                       point_colour = "black", 
                       alpha = 0.7) +
  
  # Add boxplots (smaller width to give more space to distribution)
  geom_boxplot(width = 0.12, outlier.shape = NA, alpha = 0.7, 
               fill = "white", color = "black", size = 1) +
  
  # Add a vertical line at 0 (benchmark median)
  geom_vline(aes(xintercept = 0), color = "black", size = 2) +
  
  # Custom color scale where negative values are red and positive values are blue
  scale_fill_manual(values = c("grey", "#DC0000FF") )+
                      
    # Add text annotations for the underestimation rate (place them to the right of the plot)
  geom_text(data = underestimation_rate, 
           aes(x = -50.25, y = model, label = paste0(round(under_rate, 1), "%")),  # Round percentage to 1 decimal
            vjust = -2, size = 5, aspect = 0.2, color = "#DC0000FF") +  # Adjust text size and positioning
  
  # Labels and annotation
  labs(x = "Model Error (benchmark - simulated) in % yield anomalies", y = "") +
  
  # Add facets for different climate extremes
  facet_wrap(~ climate_extreme, scales = "free_x", ncol = 2) +
  
  # Custom x-axis breaks
  scale_x_continuous(breaks = seq(-100, 100, by = 20)) +
  
  # Truncate the x-axis to focus on the range
  coord_cartesian(xlim = c(-100, 100)) +

  theme(legend.position = "none", 
        axis.title.x = element_text(size = 18),          # Increase x-axis title size
        axis.text.x = element_text(size = 18),           # Increase x-axis text size
        axis.text.y = element_text(size = 18),           # Increase y-axis text size
        strip.text = element_text(size = 18))            # Increase facet label size

raincloud_perc

Code

ggsave(file.path(base_path, "suppfig25_RaincloudExtremes.png"), raincloud_perc, width = 18, height = 20, dpi = 300)

WHEAT

Code

underestimation_rate <- crop_data[["wwh"]] %>%
mutate(error_perc = (divtrend_obs - 1)*100 - (divtrend_sim - 1)*100) %>%
  mutate(error_perc = ifelse(is.infinite(error_perc), -1000, error_perc)) %>%  # Replace Inf with -1000
  group_by(model, climate_extreme) %>%
  summarise(under_rate = mean(error_perc < 0) * 100)  # Calculate % of values below 0


# Create the raincloud plot with Inf values filtered out
raincloud_perc <- crop_data[["wwh"]] %>%
mutate(error_perc = (divtrend_obs - 1)*100 - (divtrend_sim - 1)*100) %>%
  filter(is.finite(error_perc), abs(error_perc) < 1000) %>%  # Filter out Inf values and keep within range
  mutate(model = factor(model, levels = order)) %>%  # Order the models using the defined order
  ggplot(aes(x = error_perc, y = model)) +
  coord_trans(x = scales::modulus_trans(p = 0.5)) +  # Adjust p for stronger compression
  
  # Add half-eye plot with fill based on whether the error is negative or positive
  ggdist::stat_halfeye(aes(fill = stat(x < 0)),  # Logical fill based on error sign
                       adjust = 0.7, 
                       width = 0.7,               # Increase distribution width
                       justification = -0.1,    # Move half-eye plot further from boxplot
                       .width = c(0.3, 0), 
                       point_colour = "black", 
                       alpha = 0.7) +
  
  # Add boxplots (smaller width to give more space to distribution)
  geom_boxplot(width = 0.12, outlier.shape = NA, alpha = 0.7, 
               fill = "white", color = "black", size = 1) +
  
  # Add a vertical line at 0 (benchmark median)
  geom_vline(aes(xintercept = 0), color = "black", size = 2) +
  
  # Custom color scale where negative values are red and positive values are blue
  scale_fill_manual(values = c("grey", "#DC0000FF") )+
                      
    # Add text annotations for the underestimation rate (place them to the right of the plot)
  geom_text(data = underestimation_rate, 
           aes(x = -50.25, y = model, label = paste0(round(under_rate, 1), "%")),  # Round percentage to 1 decimal
            vjust = -2, size = 5, aspect = 0.2, color = "#DC0000FF") +  # Adjust text size and positioning
  
  # Labels and annotation
  labs(x = "Model Error (benchmark - simulated) in % yield anomalies", y = "") +
  
  # Add facets for different climate extremes
  facet_wrap(~ climate_extreme, scales = "free_x", ncol = 2) +
  
  # Custom x-axis breaks
  scale_x_continuous(breaks = seq(-100, 100, by = 20)) +
  
  # Truncate the x-axis to focus on the range
  coord_cartesian(xlim = c(-100, 100)) +

  theme(legend.position = "none", 
        axis.title.x = element_text(size = 18),          # Increase x-axis title size
        axis.text.x = element_text(size = 18),           # Increase x-axis text size
        axis.text.y = element_text(size = 18),           # Increase y-axis text size
        strip.text = element_text(size = 18))            # Increase facet label size

raincloud_perc

Code

ggsave(file.path(base_path, "suppfig26_RaincloudExtremes.png"), raincloud_perc, width = 18, height = 20, dpi = 300)

SOY

Code

underestimation_rate <- crop_data[["soy"]] %>%
mutate(error_perc = (divtrend_obs - 1)*100 - (divtrend_sim - 1)*100) %>%
  mutate(error_perc = ifelse(is.infinite(error_perc), -1000, error_perc)) %>%  # Replace Inf with -1000
  group_by(model, climate_extreme) %>%
  summarise(under_rate = mean(error_perc < 0) * 100)  # Calculate % of values below 0


# Create the raincloud plot with Inf values filtered out
raincloud_perc <- crop_data[["soy"]] %>%
mutate(error_perc = (divtrend_obs - 1)*100 - (divtrend_sim - 1)*100) %>%
  filter(is.finite(error_perc), abs(error_perc) < 1000) %>%  # Filter out Inf values and keep within range
  mutate(model = factor(model, levels = order)) %>%  # Order the models using the defined order
  ggplot(aes(x = error_perc, y = model)) +
  coord_trans(x = scales::modulus_trans(p = 0.5)) +  # Adjust p for stronger compression
  
  # Add half-eye plot with fill based on whether the error is negative or positive
  ggdist::stat_halfeye(aes(fill = stat(x < 0)),  # Logical fill based on error sign
                       adjust = 0.7, 
                       width = 0.7,               # Increase distribution width
                       justification = -0.1,    # Move half-eye plot further from boxplot
                       .width = c(0.3, 0), 
                       point_colour = "black", 
                       alpha = 0.7) +
  
  # Add boxplots (smaller width to give more space to distribution)
  geom_boxplot(width = 0.12, outlier.shape = NA, alpha = 0.7, 
               fill = "white", color = "black", size = 1) +
  
  # Add a vertical line at 0 (benchmark median)
  geom_vline(aes(xintercept = 0), color = "black", size = 2) +
  
  # Custom color scale where negative values are red and positive values are blue
  scale_fill_manual(values = c("grey", "#DC0000FF") )+
                      
    # Add text annotations for the underestimation rate (place them to the right of the plot)
  geom_text(data = underestimation_rate, 
           aes(x = -50.25, y = model, label = paste0(round(under_rate, 1), "%")),  # Round percentage to 1 decimal
            vjust = -2, size = 5, aspect = 0.2, color = "#DC0000FF") +  # Adjust text size and positioning
  
  # Labels and annotation
  labs(x = "Model Error (benchmark - simulated) in % yield anomalies", y = "") +
  
  # Add facets for different climate extremes
  facet_wrap(~ climate_extreme, scales = "free_x", ncol = 2) +
  
  # Custom x-axis breaks
  scale_x_continuous(breaks = seq(-100, 100, by = 20)) +
  
  # Truncate the x-axis to focus on the range
  coord_cartesian(xlim = c(-100, 100)) +

  theme(legend.position = "none", 
        axis.title.x = element_text(size = 18),          # Increase x-axis title size
        axis.text.x = element_text(size = 18),           # Increase x-axis text size
        axis.text.y = element_text(size = 18),           # Increase y-axis text size
        strip.text = element_text(size = 18))            # Increase facet label size

raincloud_perc

Code

ggsave(file.path(base_path, "suppfig27_RaincloudExtremes.png"), raincloud_perc, width = 18, height = 20, dpi = 300)

MAIZE

Code

underestimation_rate <- crop_data[["mai"]] %>%
mutate(error_perc = (divtrend_obs - 1)*100 - (divtrend_sim - 1)*100) %>%
  mutate(error_perc = ifelse(is.infinite(error_perc), -1000, error_perc)) %>%  # Replace Inf with -1000
  group_by(model, climate_extreme) %>%
  summarise(under_rate = mean(error_perc < 0) * 100)  # Calculate % of values below 0


# Create the raincloud plot with Inf values filtered out
raincloud_perc <- crop_data[["mai"]] %>%
mutate(error_perc = (divtrend_obs - 1)*100 - (divtrend_sim - 1)*100) %>%
  filter(is.finite(error_perc), abs(error_perc) < 1000) %>%  # Filter out Inf values and keep within range
  mutate(model = factor(model, levels = order)) %>%  # Order the models using the defined order
  ggplot(aes(x = error_perc, y = model)) +
  coord_trans(x = scales::modulus_trans(p = 0.5)) +  # Adjust p for stronger compression
  
  # Add half-eye plot with fill based on whether the error is negative or positive
  ggdist::stat_halfeye(aes(fill = stat(x < 0)),  # Logical fill based on error sign
                       adjust = 0.7, 
                       width = 0.7,               # Increase distribution width
                       justification = -0.1,    # Move half-eye plot further from boxplot
                       .width = c(0.3, 0), 
                       point_colour = "black", 
                       alpha = 0.7) +
  
  # Add boxplots (smaller width to give more space to distribution)
  geom_boxplot(width = 0.12, outlier.shape = NA, alpha = 0.7, 
               fill = "white", color = "black", size = 1) +
  
  # Add a vertical line at 0 (benchmark median)
  geom_vline(aes(xintercept = 0), color = "black", size = 2) +
  
  # Custom color scale where negative values are red and positive values are blue
  scale_fill_manual(values = c("grey", "#DC0000FF") )+
                      
    # Add text annotations for the underestimation rate (place them to the right of the plot)
  geom_text(data = underestimation_rate, 
           aes(x = -50.25, y = model, label = paste0(round(under_rate, 1), "%")),  # Round percentage to 1 decimal
            vjust = -2, size = 5, aspect = 0.2, color = "#DC0000FF") +  # Adjust text size and positioning
  
  # Labels and annotation
  labs(x = "Model Error (benchmark - simulated) in % yield anomalies", y = "") +
  
  # Add facets for different climate extremes
  facet_wrap(~ climate_extreme, scales = "free_x", ncol = 2) +
  
  # Custom x-axis breaks
  scale_x_continuous(breaks = seq(-100, 100, by = 20)) +
  
  # Truncate the x-axis to focus on the range
  coord_cartesian(xlim = c(-100, 100)) +

  theme(legend.position = "none", 
        axis.title.x = element_text(size = 18),          # Increase x-axis title size
        axis.text.x = element_text(size = 18),           # Increase x-axis text size
        axis.text.y = element_text(size = 18),           # Increase y-axis text size
        strip.text = element_text(size = 18))            # Increase facet label size

raincloud_perc

Code

ggsave(file.path(base_path, "suppfig28_RaincloudExtremes.png"), raincloud_perc, width = 18, height = 20, dpi = 300)

Supplementary figures 29-33 – rainclouds for absolute model error

AGGREGATED

Code

underestimation_rate <- crop_data[["aggr"]] %>%
  mutate(error = difftrend_obs - difftrend_sim) %>%  # Replace Inf with -1000
  group_by(model, climate_extreme) %>%
  summarise(under_rate = mean(error < 0) * 100)  # Calculate % of values below 0


# Create the raincloud plot with Inf values filtered out
raincloud_perc <- crop_data[["aggr"]] %>%
mutate(error = difftrend_obs - difftrend_sim) %>%
  mutate(model = factor(model, levels = order)) %>%  # Order the models using the defined order
  ggplot(aes(x = error, y = model)) +
  coord_trans(x = scales::modulus_trans(p = 0.5)) +  # Adjust p for stronger compression
  
  # Add half-eye plot with fill based on whether the error is negative or positive
  ggdist::stat_halfeye(aes(fill = stat(x < 0)),  # Logical fill based on error sign
                       adjust = 0.7, 
                       width = 0.7,               # Increase distribution width
                       justification = -0.1,    # Move half-eye plot further from boxplot
                       .width = c(0.3, 0), 
                       point_colour = "black", 
                       alpha = 0.7) +
  
  # Add boxplots (smaller width to give more space to distribution)
  geom_boxplot(width = 0.12, outlier.shape = NA, alpha = 0.7, 
               fill = "white", color = "black", size = 1) +
  
  # Add a vertical line at 0 (benchmark median)
  geom_vline(aes(xintercept = 0), color = "black", size = 2) +
  
  # Custom color scale where negative values are red and positive values are blue
  scale_fill_manual(values = c("grey", "#DC0000FF") )+
                      
    # Add text annotations for the underestimation rate (place them to the right of the plot)
  geom_text(data = underestimation_rate, 
           aes(x = -1.25, y = model, label = paste0(round(under_rate, 1), "%")),  # Round percentage to 1 decimal
            vjust = -2, size = 5, aspect = 0.2, color = "#DC0000FF") +  # Adjust text size and positioning
  
  # Labels and annotation
  labs(x = "Model Error (benchmark - simulated) in tons/ha yield anomalies", y = "") +
  
  # Add facets for different climate extremes
  facet_wrap(~ climate_extreme, scales = "free_x", ncol = 2) +
  
  # Custom x-axis breaks
  scale_x_continuous(breaks = seq(-3, 3, by = 1)) +
  
  # Truncate the x-axis to focus on the range
  coord_cartesian(xlim = c(-3, 3)) +

  theme(legend.position = "none", 
        axis.title.x = element_text(size = 18),          # Increase x-axis title size
        axis.text.x = element_text(size = 18),           # Increase x-axis text size
        axis.text.y = element_text(size = 18),           # Increase y-axis text size
        strip.text = element_text(size = 18))            # Increase facet label size

raincloud_perc

Code

ggsave(file.path(base_path, "suppfig29_RaincloudExtremes.png"), raincloud_perc, width = 18, height = 20, dpi = 300)

RICE

Code

underestimation_rate <- crop_data[["ri1"]] %>%
  mutate(error = difftrend_obs - difftrend_sim) %>%  # Replace Inf with -1000
  group_by(model, climate_extreme) %>%
  summarise(under_rate = mean(error < 0) * 100)  # Calculate % of values below 0


# Create the raincloud plot with Inf values filtered out
raincloud_perc <- crop_data[["ri1"]] %>%
mutate(error = difftrend_obs - difftrend_sim) %>%
  mutate(model = factor(model, levels = order)) %>%  # Order the models using the defined order
  ggplot(aes(x = error, y = model)) +
  coord_trans(x = scales::modulus_trans(p = 0.5)) +  # Adjust p for stronger compression
  
  # Add half-eye plot with fill based on whether the error is negative or positive
  ggdist::stat_halfeye(aes(fill = stat(x < 0)),  # Logical fill based on error sign
                       adjust = 0.7, 
                       width = 0.7,               # Increase distribution width
                       justification = -0.1,    # Move half-eye plot further from boxplot
                       .width = c(0.3, 0), 
                       point_colour = "black", 
                       alpha = 0.7) +
  
  # Add boxplots (smaller width to give more space to distribution)
  geom_boxplot(width = 0.12, outlier.shape = NA, alpha = 0.7, 
               fill = "white", color = "black", size = 1) +
  
  # Add a vertical line at 0 (benchmark median)
  geom_vline(aes(xintercept = 0), color = "black", size = 2) +
  
  # Custom color scale where negative values are red and positive values are blue
  scale_fill_manual(values = c("grey", "#DC0000FF") )+
                      
    # Add text annotations for the underestimation rate (place them to the right of the plot)
  geom_text(data = underestimation_rate, 
           aes(x = -1.25, y = model, label = paste0(round(under_rate, 1), "%")),  # Round percentage to 1 decimal
            vjust = -2, size = 5, aspect = 0.2, color = "#DC0000FF") +  # Adjust text size and positioning
  
  # Labels and annotation
  labs(x = "Model Error (benchmark - simulated) in tons/ha yield anomalies", y = "") +
  
  # Add facets for different climate extremes
  facet_wrap(~ climate_extreme, scales = "free_x", ncol = 2) +
  
  # Custom x-axis breaks
  scale_x_continuous(breaks = seq(-3, 3, by = 1)) +
  
  # Truncate the x-axis to focus on the range
  coord_cartesian(xlim = c(-3, 3)) +

  theme(legend.position = "none", 
        axis.title.x = element_text(size = 18),          # Increase x-axis title size
        axis.text.x = element_text(size = 18),           # Increase x-axis text size
        axis.text.y = element_text(size = 18),           # Increase y-axis text size
        strip.text = element_text(size = 18))            # Increase facet label size

raincloud_perc

Code

ggsave(file.path(base_path, "suppfig30_RaincloudExtremes.png"), raincloud_perc, width = 18, height = 20, dpi = 300)

WHEAT

Code

underestimation_rate <- crop_data[["wwh"]] %>%
  mutate(error = difftrend_obs - difftrend_sim) %>%  # Replace Inf with -1000
  group_by(model, climate_extreme) %>%
  summarise(under_rate = mean(error < 0) * 100)  # Calculate % of values below 0


# Create the raincloud plot with Inf values filtered out
raincloud_perc <- crop_data[["wwh"]] %>%
mutate(error = difftrend_obs - difftrend_sim) %>%
  mutate(model = factor(model, levels = order)) %>%  # Order the models using the defined order
  ggplot(aes(x = error, y = model)) +
  coord_trans(x = scales::modulus_trans(p = 0.5)) +  # Adjust p for stronger compression
  
  # Add half-eye plot with fill based on whether the error is negative or positive
  ggdist::stat_halfeye(aes(fill = stat(x < 0)),  # Logical fill based on error sign
                       adjust = 0.7, 
                       width = 0.7,               # Increase distribution width
                       justification = -0.1,    # Move half-eye plot further from boxplot
                       .width = c(0.3, 0), 
                       point_colour = "black", 
                       alpha = 0.7) +
  
  # Add boxplots (smaller width to give more space to distribution)
  geom_boxplot(width = 0.12, outlier.shape = NA, alpha = 0.7, 
               fill = "white", color = "black", size = 1) +
  
  # Add a vertical line at 0 (benchmark median)
  geom_vline(aes(xintercept = 0), color = "black", size = 2) +
  
  # Custom color scale where negative values are red and positive values are blue
  scale_fill_manual(values = c("grey", "#DC0000FF") )+
                      
    # Add text annotations for the underestimation rate (place them to the right of the plot)
  geom_text(data = underestimation_rate, 
           aes(x = -1.25, y = model, label = paste0(round(under_rate, 1), "%")),  # Round percentage to 1 decimal
            vjust = -2, size = 5, aspect = 0.2, color = "#DC0000FF") +  # Adjust text size and positioning
  
  # Labels and annotation
  labs(x = "Model Error (benchmark - simulated) in tons/ha yield anomalies", y = "") +
  
  # Add facets for different climate extremes
  facet_wrap(~ climate_extreme, scales = "free_x", ncol = 2) +
  
  # Custom x-axis breaks
  scale_x_continuous(breaks = seq(-3, 3, by = 1)) +
  
  # Truncate the x-axis to focus on the range
  coord_cartesian(xlim = c(-3, 3)) +

  theme(legend.position = "none", 
        axis.title.x = element_text(size = 18),          # Increase x-axis title size
        axis.text.x = element_text(size = 18),           # Increase x-axis text size
        axis.text.y = element_text(size = 18),           # Increase y-axis text size
        strip.text = element_text(size = 18))            # Increase facet label size

raincloud_perc

Code

ggsave(file.path(base_path, "suppfig31_RaincloudExtremes.png"), raincloud_perc, width = 18, height = 20, dpi = 300)

SOY

Code

underestimation_rate <- crop_data[["soy"]] %>%
  mutate(error = difftrend_obs - difftrend_sim) %>%  # Replace Inf with -1000
  group_by(model, climate_extreme) %>%
  summarise(under_rate = mean(error < 0) * 100)  # Calculate % of values below 0


# Create the raincloud plot with Inf values filtered out
raincloud_perc <- crop_data[["soy"]] %>%
mutate(error = difftrend_obs - difftrend_sim) %>%
  mutate(model = factor(model, levels = order)) %>%  # Order the models using the defined order
  ggplot(aes(x = error, y = model)) +
  coord_trans(x = scales::modulus_trans(p = 0.5)) +  # Adjust p for stronger compression
  
  # Add half-eye plot with fill based on whether the error is negative or positive
  ggdist::stat_halfeye(aes(fill = stat(x < 0)),  # Logical fill based on error sign
                       adjust = 0.7, 
                       width = 0.7,               # Increase distribution width
                       justification = -0.1,    # Move half-eye plot further from boxplot
                       .width = c(0.3, 0), 
                       point_colour = "black", 
                       alpha = 0.7) +
  
  # Add boxplots (smaller width to give more space to distribution)
  geom_boxplot(width = 0.12, outlier.shape = NA, alpha = 0.7, 
               fill = "white", color = "black", size = 1) +
  
  # Add a vertical line at 0 (benchmark median)
  geom_vline(aes(xintercept = 0), color = "black", size = 2) +
  
  # Custom color scale where negative values are red and positive values are blue
  scale_fill_manual(values = c("grey", "#DC0000FF") )+
                      
    # Add text annotations for the underestimation rate (place them to the right of the plot)
  geom_text(data = underestimation_rate, 
           aes(x = -1.25, y = model, label = paste0(round(under_rate, 1), "%")),  # Round percentage to 1 decimal
            vjust = -2, size = 5, aspect = 0.2, color = "#DC0000FF") +  # Adjust text size and positioning
  
  # Labels and annotation
  labs(x = "Model Error (benchmark - simulated) in tons/ha yield anomalies", y = "") +
  
  # Add facets for different climate extremes
  facet_wrap(~ climate_extreme, scales = "free_x", ncol = 2) +
  
  # Custom x-axis breaks
  scale_x_continuous(breaks = seq(-3, 3, by = 1)) +
  
  # Truncate the x-axis to focus on the range
  coord_cartesian(xlim = c(-3, 3)) +

  theme(legend.position = "none", 
        axis.title.x = element_text(size = 18),          # Increase x-axis title size
        axis.text.x = element_text(size = 18),           # Increase x-axis text size
        axis.text.y = element_text(size = 18),           # Increase y-axis text size
        strip.text = element_text(size = 18))            # Increase facet label size

raincloud_perc

Code

ggsave(file.path(base_path, "suppfig32_RaincloudExtremes.png"), raincloud_perc, width = 18, height = 20, dpi = 300)

MAIZE

Code

underestimation_rate <- crop_data[["mai"]] %>%
  mutate(error = difftrend_obs - difftrend_sim) %>%  # Replace Inf with -1000
  group_by(model, climate_extreme) %>%
  summarise(under_rate = mean(error < 0) * 100)  # Calculate % of values below 0


# Create the raincloud plot with Inf values filtered out
raincloud_perc <- crop_data[["mai"]] %>%
mutate(error = difftrend_obs - difftrend_sim) %>%
  mutate(model = factor(model, levels = order)) %>%  # Order the models using the defined order
  ggplot(aes(x = error, y = model)) +
  coord_trans(x = scales::modulus_trans(p = 0.5)) +  # Adjust p for stronger compression
  
  # Add half-eye plot with fill based on whether the error is negative or positive
  ggdist::stat_halfeye(aes(fill = stat(x < 0)),  # Logical fill based on error sign
                       adjust = 0.7, 
                       width = 0.7,               # Increase distribution width
                       justification = -0.1,    # Move half-eye plot further from boxplot
                       .width = c(0.3, 0), 
                       point_colour = "black", 
                       alpha = 0.7) +
  
  # Add boxplots (smaller width to give more space to distribution)
  geom_boxplot(width = 0.12, outlier.shape = NA, alpha = 0.7, 
               fill = "white", color = "black", size = 1) +
  
  # Add a vertical line at 0 (benchmark median)
  geom_vline(aes(xintercept = 0), color = "black", size = 2) +
  
  # Custom color scale where negative values are red and positive values are blue
  scale_fill_manual(values = c("grey", "#DC0000FF") )+
                      
    # Add text annotations for the underestimation rate (place them to the right of the plot)
  geom_text(data = underestimation_rate, 
           aes(x = -1.25, y = model, label = paste0(round(under_rate, 1), "%")),  # Round percentage to 1 decimal
            vjust = -2, size = 5, aspect = 0.2, color = "#DC0000FF") +  # Adjust text size and positioning
  
  # Labels and annotation
  labs(x = "Model Error (benchmark - simulated) in tons/ha yield anomalies", y = "") +
  
  # Add facets for different climate extremes
  facet_wrap(~ climate_extreme, scales = "free_x", ncol = 2) +
  
  # Custom x-axis breaks
  scale_x_continuous(breaks = seq(-5, 5, by = 1)) +
  
  # Truncate the x-axis to focus on the range
  coord_cartesian(xlim = c(-5, 5)) +

  theme(legend.position = "none", 
        axis.title.x = element_text(size = 18),          # Increase x-axis title size
        axis.text.x = element_text(size = 18),           # Increase x-axis text size
        axis.text.y = element_text(size = 18),           # Increase y-axis text size
        strip.text = element_text(size = 18))            # Increase facet label size

raincloud_perc

Code

ggsave(file.path(base_path, "suppfig33_RaincloudExtremes.png"), raincloud_perc, width = 18, height = 20, dpi = 300)